Synchronization Transition in the Kuramoto Model with Colored Noise 



Ralf Tonjes 

Ochadai Academic Production, Ochanomizu University, Tokyo 112-8610, Japan 

We present a linear stability analysis of the incoherent state in a system of globally coupled, 
identical phase oscillators subject to colored noise. In that we succeed to bridge the extreme time 
scales between the formerly studied and analytically solvable cases of white noise and quenched 
random frequencies. 
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The term Kuramoto Model refers to a class of nonlin- 
ear models which describe the dynamics of autonomous 
limit cycle oscillators by phase equations and interac- 
tions between them via coupling functions of phase dif- 
ferences. Since its original formulation [l| it has been 
modified to include for instance other nonlinear effects, 
coupling topology or delayed coupling An important 
property of these models is the existence of a transition 
to synchronization in large systems of coupled oscillators 
mediated through the opposing effects of attractive inter- 
action and heterogeneity. Synchronization is a collective 
phenomenon where the phases of the oscillators become 
correlated leading to macroscopic oscillations or more 
complicated behavior ^S-J]- The Kuramoto Model is 
therefore able to reproduce a fundamental mechanism of 
self-organization in nature which is important for pattern 
formation, information processing and transport among 
others. 

In Kuramoto considers the case of all-to-all coupling 
where each oscillator couples equally strong to all other 
oscillators in the system. The Kuramoto phase equations 
for such a system are 
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where 'dn is the phase of the oscillator n, arjn is an indi- 
vidual force which may be the natural frequency of the 
oscillator or a time dependent perturbation, g{Ad) is a 
periodic coupling function of a phase difference and N 
is the total number of oscillators. Disorder is realized 
through a distribution of random forces crrjn, where a de- 
notes the noise amplitude in units of coupling strength. 
When the forces are time independent the system mod- 
els an ensemble of oscillators with nonidentical natural 
frequencies. For quenched random frequencies with uni- 
modal distribution a continuous phase transition from 
an incoherent regime of evenly distributed phases to a 
regime of partial synchronization can be observed when 
a is changed. In fact, depending on the shape of the fre- 
quency distribution or the couplingfunction, even more 
complicated behavior is possible 0,11, |1] . If, on the other 
hand, the r]n change very rapidly, they may be approxi- 
mated by white noise. Again, a continuous phase tran- 
sition is predicted as the noise strength is changed Q. 



These two analytically solvable cases mark the extremes 
of time scale separation between the dynamics of the os- 
cillators and the fluctuations. In experiments, however, 
system parameters may drift at time scales comparable 
to the drift of the oscillator phase differences. Moreover, 
if the random forces arjn are intrinsic to the system, for 
instance, in phase coherent chaotic oscillators, or in a 
random network of identical phase oscillators 0, , the 
time scales are not necessarily separated. It is there- 
fore of great interest to know how the critical coupling 
strength for the phase transition to partial synchroniza- 
tion for colored noise differs from the known values at 
quenched or white noise. 

Numerical investigations of that question have been car- 
ried out and qualitative answers have been given at se- 
lected parameters 
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However, an analytical solution 
to the problem has so far remained an open problem. 
Here we provide the solution in the two cases of the ran- 
dom telegraph process (TP) and the Ornstein-Uhlenbeck 
process (OU) as source of the colored noise. We find, that 
the type of the random process is essential for the tran- 
sition point to synchronization, as can be expected from 
the rich behavior of the Kuramoto model with different 
quenched frequency distributions |3, [13] . 



Evolution of Phase and Frequency Distribution 

In the thermodynamical limit N ~> oo the system can 
be described by a density r], t) of phases -d and forces 
rj. The evolution of this density is given by a Fokker- 
Planck equation 



dtp = L{, [p] p -f L^p 



(2) 



where we assume that the forces rjn are independent, lin- 
ear random processes described by a linear operator L^, 
whereas the Fokker-Planck operator [p] ■ = —d^{'d-) 
that acts on the phases depends on the mean field and 
is therefore a functional of the oscillator density p. Our 
strategy is to linearize Eq. ^ around the stationary in- 
coherent state p{'d, rj, t) = p{t])/2tt and look for a critical 
condition of the stability of the eigenmodes of p. If has 
a finite number of eigenmodes, as is the case for the ran- 
dom telegraph process, we will only have to solve a finite 
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system of linear equations. This is not the case, how- 
ever, when we consider the probably most applied case 
of the OU process. Then has an infinite but countable 
number of eigenmodes and we are faced to determine the 
stability of an infinite system of linear ODEs. 
Given the eigenvalues A„ and eigenmodes ifiniv) of -^») 
with Lfjipn — Xn^n and = Aq > ReAi ... we start by 
expanding p{'d, if) and g(A-!5) as 

fc= — oo n=0 

(3) 
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The Kuramoto phase equation ([1]) gives 

/TT nOO 
'TT J — OO 

(4) 

OO 

= <jrj+ J2 ^fco g*k e'" . 

fc — — OO 

Inserting Eqs. ([H} and ([4]) into we find for the modes 
Zkn the nonlinear equation 

OO OO 

Zkn = ^ ZiogiZ(^k-l)7i - '^kcr'^ ZkmMm7i-\-Kzkn 

1— — 00 m—0 

(5) 

with 

OO 

T]ip„{r]) = ^ MnrnVmiv) ■ (6) 
m—0 

We remark that the celebrated ansatz of Antonsen and 
Ott [1] corresponding to Ckiv) '■= J2n ^kn^fin / ^fio = a^{ri) 
(a*l'^l for k <Q) only leads to a closed nonlinear dynamic 
equation for a(ry) when all eigenvalues An are equal to 
zero, i.e., in the limit of quenched random frequencies. 
The order parameter R = | (e*'') | is equal to the abso- 
lute value of zio which is zero when the phases are dis- 
tributed uniformly in the incoherent state. One can check 
that the incoherent distribution 77, = (pQ{r])/2Tr or 
Zkn = ^fco<^nO IS a Stationary solution of the Fokker-Planck 
equation Thus, keeping only the terms linear in the 
small quantities z^n we obtain linearized equations for 
the dynamics of the modes 

Zko = -ifc (ffo -I- gl) Zko - ika ZkmMmo, 

m 

(7) 

Zkn = (A„ - ifc^o) ^kn - '^k(J ^ ZkmMmn- 

m 

The Fourier modes of the probability distribution, i.e. 
the Zkn for different fc, do not interact linearly with 
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FIG. 1: (color online) Case of the random telegraph process. 
Subfigures (a) and (b) show the parameter regions of linear 
stability of the incoherent state for (a) switching rate 7 and 
amplitude or (b) inverse amplitude and noise strength 
D-^. Shown are the two cases a — (bold dark/blue) and 
a — Tr/4 (bold light/green). The white noise limit is reached 
on the vertical axis of (b) where and — ^/2a^ = 

const. Subfigure (b) can directly be compared to the case of 
Ornstein-Uhlenbeck type noise in Fig. [21 The dashed lines are 
asymptotes obtained from Eq. ((9|. Subfigure (c) shows the 
time-averaged order parameter. The theoretical predictions 
are verified by direct numerical simulation in systems of size 
A*' = 5000 at 7 = 0.35 using as bifurcation parameter. 
The corresponding one dimensional curve is shown as dashed 
and dotted (red) line in (a) and (b). The vertical dotted lines 
mark the theoretical transition points. Hysteresis is observed 
for a = 0.0 (blue circles and squares). Subfigure (d) shows the 
phases of 1000 oscillators in the partially synchronized state 
for a — 0.0, 7 = 0.35 and = 0.3 with a phase histogram in 
the inset. In this state the oscillators are phase locked to the 
stationary mean field. 

one another and can be studied separately. The term 
— ifc^Q — — ifcgo which appears as a self-interaction term 
for all eigenmodes can be neglected for the stability anal- 
ysis since it is imaginary and only results in a bias to all 
frequencies. The incoherent state becomes unstable when 
the largest real part of the eigenvalues of the linear ODE 
([7]) becomes positive. For a linear random process with 
a finite number of eigenmodes ^Pniv) '^6 only need to to 
determine the eigenvalues of a square matrix depending 
on the system parameters. 



Random Telegraph Process 

Consider the Kuramoto model with sinusoidal coupling 
function g{Ad) — sm{A-d — a), attracting (cos a > 0) 
and with independent random forces ry„ G {—1,1} which 
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change sign as a dichotomous random Markov process 
with equal transition rate 7 between both values. Fol- 
lowing the analysis in the previous section we find that 
the linear stability of the first Fourier mode at the in- 
coherent state is determined by the eigenvalues of the 
matrix 



1(7 

-27 



(8) 



Given a and the flipping rate 7, necessary and sufficient 
conditions for stability are 



cr > 7 cos a 1 



■ 2 

sm a 



(47 — cosa)^ 



cos a -47 < 0. (9) 



The frequency fl of the mean field at the bifurcation is 
n — 27 sin a/ (47 — cos a). Interestingly increasing a the 
incoherent state may actually become unstable in some 
regions of parameter space (see Fig. [l^,b) even though it 
is stable everywhere when a — ^ 7r/2. To test our result 
we integrated Eq. ([T]) numerically with time steps de- 
termined by the random switching events and find it in 
good agreement with the prediction (Fig. [TJ;). Both the 
thermodynamic and the white noise limit are not acces- 
sible through this integration scheme since the time step 
is of order 0{1/Nj). However, the white noise limit with 
D = 2(t2/7 is recovered from Eq. ([9]) with Z?cr = 0.5 cos a 
and — 0.5 sin a. From the linear stability analysis we 
are not able to predict whether the Hopf-Bifurcation is 
supercritical or subcritical. In fact, the simulations show 
either behavior in different parameter regions. 



Ornstein-Uhlenbeck Process 

Instead of randomly switching between two values, we 
now consider i.i.d. random forces diffusing in the fashion 
of an OU process with Langevin equation 



77 = -777 + 



(10) 



and white noise {^{t)^{t')) ^ 6{t - t'). The rate 7 de- 
termines the time scale of the diffusion. To visualize 
both the white noise and the quenched noise limit it is 
of advantage to use the parameter D = cr^/7. Then the 
white noise limit with finite noise strength D is reached 
as — )■ and quenched noise corresponds to — 0. 
The eigenvalues and eigenfunctions of the Fokker-Planck 
operater for an OU process are intimately related to 
those of the quantum harmonic oscillator |13| . One finds 
A,i = -7n and M„„ = Smn-i\/n + Smn+Wn + 1- This 
turns the eigenvalue problem of Eq. ([7]) into an infinite 
system of second order difference equations. 
While determining the eigenvalues of the ODE ^ de- 
pending on D, a and a presents a major difficulty, this 
is actually not necessary. Instead we notice that at the 



transition to synchronization there is an imaginary eigen- 
value iri, which gives us an implicit condition for the bi- 
furcation line of the first Fourier mode (fc = 1) 

irizio = -igjzio - iCTZii, 

(11) 

ir2zi„ = —a^D^^nzYn — ia\fnzin~\ ^ i^V n + Izin+i- 

At the transition D, is the frequency of the emerging mean 
field. Denoting /i„ — ~i\Jn + lzi„^i/zi„ these differ- 
ence equations define a continued fraction 



/^o(a;,y) = -Kr=i 



(12) 



nx + ly 



where the dimensionless quantities x — ja~^ = D~^a 
and y = fla~^ relate the dynamical time scales in the 
system. This equation for the critical condition is one 
of the main results of this paper. The complex function 
fiQ{x,y) can be calculated efficiently from Eq. (IT^ . Us- 
ing a technique of Euler 14j , one can find the continued 
fraction in terms of functions related to confluent hyper- 
geometric functions of the first kind 



fJ'o(x,y) 



1 ifi{2,x-'{iy + x-')+2,x-^) 
X ih{l,x-^iy + x-^) + l,x-^) ' 



(13) 



1/1 (a, 6, z) = / duu''-\l-u)''-''-^e^'' 
Jo 



With igl = — cxp(ia)/2 it follows that 

(T^^ = 2|/xo(x, ?/) — iy| , and cos a = — 2(TRe/io- (14) 

The critical lines in Fig. [5] are parametrized by the time 
scale ratio x = ... 00. For fixed nonzero a the time 
scale ratio y{x) has to be determined numerically from 
Eq. p4)) . For vanishing a, /io is real, i.e. y = 0. The 
white noise limit Da- = 0.5 cos a and 57 — 0.5 sin a can 
easily be obtained from Eq. letting a —i' 00. The 
quenched noise limit a; is not trivial. For a = 
it must be compared to a^,} = 2/7ryjo(0) = a/S/tt 0). 
No simple expression exists for a =^ [l^j. As a spe- 
cial case of Eq. ()12|) . for a = and cc = 1, we obtain 
i^„.i-a-.i-2/(e-l) ;i4]. 

To test our analytic result for the critical condition 
Monte-Carlo simulations of the Kuramoto model Eq. ([1]) 
with OU random forces have been carried out with fi- 
nite step size dt. The displacements of phase 1) and force 
rj can be drawn directly from the transition probability 
p(i9 + d-d, rj + drj, t + dt\d, rj, t) under the assumption of a 
slowly changing mean field force /(??) which is assumed 
constant during an integration step. The two random 
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FIG. 2: (color online) Synchronization transition of the Ku- 
ramoto model Eq. ((T} subject to random forces of Ornstein- 
Uhlenbeck type with variance and dissipation rate 7. 
Shown is the time averaged order parameter J? as a func- 
tion of and = 7/0-'^ for (a) a — 0.0 and (b) a = 0.7 
in a system of size A'^ = 5000 averaged over 150 time units 
after a transient of 50 units. The solid white line marks the 
critical conditions Eq. (|14|l obtained by changing x — D~^a 
from 0.01 to 20. The white noise limit is located on the or- 
dinate axis for cr^^ and the quenched noise limit on the 



abscissa axis for D 



0. 



variables 



(15) 



r2 = 77(t + dt)-e-'^'^"''^*77(t) 



are Gaussian [l^l with correlation matrix 

Sn = Z?V-2(2a2z?-idi-3 + 4e-'^'^"*-e-2-'^"''^*), 

S12 - S21 = Da-\l-e-^^''''"''')\ (16) 



S22 = 1 - e 



Both values in Eq. (jTS]) can be sampled at all time 
scales and in particular also for a^^ — )■ as well as 
D^^ — >■ 0. The critical line obtained from Eq. (IT^ and 
the numerical simulations agree very well (Fig. [5]). 



By means of linear stability analysis we have suc- 
ceeded to find critical conditions for the transition to 
synchronization in the Kuramoto model of globally 
coupled, identical oscillators subject to independent but 
identically distributed colored noise forces in the cases of 
Ornstein-Uhlenbeck type and random Telegraph noise. 
We are hopeful that our results can be applied to obtain 
qualitative and quantitative predictions for the critical 
coupling strength in an ensemble of phase coherent 
chaotic oscillators 15|, Jig or networks of identical 
autonomous oscillators [7|. For such an application it 
will be necessary to approximate the experimentally 
accessible fluctuations in single oscillators by a linear 
model such as the Ornstein-Uhlenbeck process or a finite 
state Markov model like the random telegraph process. 
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